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Abstract 

Latin hypercube sampling (LHS) is generalized in terms of a spectrum of stratified sampling (SS) 
designs referred to as partially stratified sample (PSS) designs. True SS and LHS are shown to 
represent the extremes of the PSS spectrum. The variance of PSS estimates is derived along with 
some asymptotic properties. PSS designs are shown to reduce variance associated with variable 
interactions, whereas LHS reduces variance associated with main effects. Challenges associated 
with the use of PSS designs and their limitations are discussed. To overcome these challenges, 
the PSS method is coupled with a new method called Latinized stratified sampling (LSS) that 
produces sample sets that are simultaneously SS and LHS. The LSS method is equivalent to an 
Orthogonal Array based LHS under certain conditions but is easier to obtain. Utilizing an LSS on 
the subspaces of a PSS provides a sampling strategy that reduces variance associated with both 
main effects and variable interactions and can be designed specially to minimize variance for a given 
problem. Several high-dimensional numerical examples highlight the strengths and limitations of 
the method. The Latinized partially stratified sampling method is then applied to identify the best 
sample strategy for uncertainty quantification on a plate buckling problem. 

Keywords: Uncertainty Quantification, Monte Carlo Simulation, Stratified Sampling, Latin 
Hypercube Sampling 


1. Introduction 

Latin hypercube sampling (LHS) [U [2] is the most widely used random sampling method for 
Monte Carlo-based uncertainty quantification, employed in nearly every field of computational 
science, engineering, and mathematics. The seminal work by McKay et al. [T] introducing Latin 
hypercube sampling is a classic in the field of design of computer experiments. LHS is an especially 
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powerful and useful sampling method thanks primarily to the properties identified by Stein [3] who 
showed that LHS has the effect of filtering the variance associated with the additive components of 
a transformation (or main effects). This result, combined with the Hierarchical Ordering Principle 
[3] - which states that main effects and low order interactions are likely to govern most general 
transformations - causes LHS to reduce variance significantly for many applications. 

The widespread popularity of LHS has led to the invention of numerous variants meant to im¬ 
prove space-filling [siiiiTiisiEKinKn], optimize projective properties [12], minimize least square 
error and maximize entropy |13j . and reduce spurious correlations [HKiaiisiiiTiiTiiiiiiHi. Mean¬ 
while, LHS has been applied to nearly every type of probabilistic analysis one can imagine, ranging 
from the estimation of reliability (probability of failure) [T9| to coefficient estimation for polynomial 
chaos, neural network, and other types of surrogate models [Ml [ZD. The intention of this paper 
is not to present another variant of the LHS methodology or to apply it in a new or novel way. 
Rather, we present a broad generalization of the methodology in the context of stratified sampling 
- from which LHS is derived. 

Stratified sampling (SS), the “parent” methodology of LHS has been widely used in the social 
sciences and financial mathematics owing to its ability to partition a population into strata (or 
categories) that can be weighted according to their conditional probabilities. Some recent develop¬ 
ments have begun to encourage its use in Monte Carlo analysis of computer models by adaptively 
stratifying the probability space of a random vector |22l [Z5] . These methods rely on “true” strati¬ 
fied sampling wherein all dimensions of the space are stratified simultaneously allowing the analyst 
to concentrate samples in probabilistically weighted regions of the space that are important for the 
problem at hand. LHS, meanwhile, lies at the opposite end of the “spectrum” of stratified sampling 
methods (Figure such that each dimension of the random vector is stratified individually and 
the vector is constructed through random pairings. 

In this work, the intermediate space on the spectrum of stratified sampling methods is explored 
such that stratihcation can occur on any set of Aij-dimensional orthogonal subspaces of the N- 
dimensional sample space S subject to Ni = N. This generalized sampling method is referred 
to as Partially Stratihed Sampling (PSS) and the properties of their designs are explored in detail. 
In particular, the PSS method is shown to reduce variance associated with low-dimensional inter¬ 
actions within a high-dimensional transformation. Some discussion surrounding the appropriate 
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Figure 1: Spectrum of stratified sampling methods. 


use of a PSS design is provided and a powerful hybrid PSS-LHS method is proposed that simulta¬ 
neously reduces variance associated with the main-effects and low-dimensional interactions. This 
method, referred to as “Latinized” Partially Stratihed Sampling (LPSS), combines the variance 
reductions of both PSS and LHS to yield a major improvement in sample efficiency. Several high¬ 
dimensional demonstration problems are presented and the method is applied to the probabilistic 
assessment of plate buckling strength where the interaction of geometric and material variables are 
very important. 

2. Review^ of Sampling Methods 

This section provides a brief review of the sampling methods used for analysis in this work. 
We will consider only uncorrelated random variables as it is common practice to maps correlated 
variables onto a set of uncorrelated ones using, for example. Principal Component Analysis. 

2.1. Simple Random Sampling 

Classical Monte Carlo methods rely on so-called Simple Random Sampling (SRS) or Monte Carlo 
Sampling in which realizations of the vector x (samples) are generated as independent and iden¬ 
tically distributed (iid) realizations on S with marginal cumulative distribution functions (CDFs) 
Dxii-) by: 

Xi = D^^^{Ui)-,i = 1,2,...,n ( 1 ) 

where Ut are iid uniformly distributed samples on [0,1]. The realizations x are then applied to the 
system y = F(x) and y is statistically evaluated. 
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2.2. Stratified Sampling 

Stratified Sampling begins by dividing the sample space S into a collection of M disjoint subsets 
(strata) = 1,2, ...,M with = S and fip n Q,q = %]p q. Sample realizations 

from a given stratum k, x^. = {xik,X 2 k, ■ ■ ■ ,X]\fk}, are generated by randomly sampling the vector 
components according to: 

Xik = ■ ,N (2) 

where Uik are iid uniformly distributed samples on with = DxfiClk) ^ik = ^XiiQk)^ 

and Cifc and denote the lower and upper bounds respectively of the vector component of 
stratum Typically the stratification is performed directly in the probability space meaning 
that the strata are defined directly by the bounds and 

2.3. Latin Hypercube Sampling 

Latin Hypercube Sampling (LHS) operates by dividing the subspace of each vector component 
5*; i = 1, 2,..., into M = n disjoint subsets (strata) of equal probability Ltik] i = 1,2,..., N; k = 
1,2,..., M. Samples of each vector component are drawn from the respective strata according to: 

Xik — ^Xi i^ik): i — 1,2,..., N ',k = 1,2,..., M (3) 

where Uik are iid uniformly distributed samples on with ^\. = ^ and ^k ~ 

samples x are assembled by randomly grouping the terms of the generated vector components. 
That is, a term Xik is randomly selected from each vector component (without replacement) and 
these terms are grouped to produce a sample. This process is repeated M = n times. 

Because the component samples are randomly paired, an LHS is not unique; there are 
possible combinations. With this in mind, improved LHS algorithms iterate to determine optimal 
pairings according to some specihed criteria - such as reduced correlation among the terms or 
enhanced space-filling properties (e.g. ismsidziiiainiEiiin]). 

2 . 4 . Variance reduction in stratified sampling and Latin hypercube sampling 
Consider the general statistical estimator defined by: 

n 

T{yi,...,yn) = '^wig{yi) (4) 

i=i 

where yi = h{xi) and xi denotes a sample generated using SRS, SS, or LHS, wi are sample weights, 
and g{-) is an arbitrary function. Note that when g{y) = y'^, T is an estimate of the r^^ moment 
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and when g(y) = I{y < Y}, where /{•} denotes the indicator function, T is the empirical CDF. For 
conventional Monte Carlo analysis, wi = -^ \/l and the variance of the statistical estimator {Tr) is 
given by Var[rij] = ^ where =Var[gf(y)]. Classical Monte Carlo estimates generally serve as 
the measure by which variance reduction techniques are compared. 

SS and LHS are both techniques to reduce the variance of statistical estimators when compared 
to classical Monte Carlo estimates although they do so through different statistical mechanisms. 
Stratified sampling has been proven to unconditionally reduce the variance of statistical estimators 
(denoted T 5 ) when compared to SRS such that the variance reduction depends on the differences 
between the strata means gk and the overall mean r as [I]: 

1 ^ 

Var[T 5 ] = VarfFij]- '^PkifJ-k - rf (5) 

^ k=l 

when the strata are sampled proportionately such that the number of samples on stratum k, 
nk = Pk^- LHS, on the other hand, reduces variance by creating negative covariance between 
sample cells such that the variance of a LHS estimator can be expressed as [T]: 

Var[rL] = Varfr^] + 

R 

where r is, again, the mean value, pi is the mean of LHS cell i, and R denotes the n^{n — 1)'^ pairs 
{pi, pj) of cells having no cell coordinates in common. Note that the variance is reduced only when 
the second term (cell covariance) is negative. This, however, is not a restrictive condition as McKay 
et al. [T] have shown that variance is always reduced for monotonic functions Y = h(X.) and Stein 
[3] showed further that variance is reduced for any function /i(-) having finite second moment when 
n >> N. Moreover, Stein showed that LHS reduces variance by effectively filtering the main effects 
of the transformation such that the stronger the main effects, the more the variance is reduced. 
This will be discussed in more detail later. 

3. Partially Stratified Sampling 

True stratified sampling and Latin hypercube sampling can be viewed as two extremes on a 
spectrum of possible stratification methods (Figure]^. In true SS, the stratification occurs on 
all dimensions simultaneously while LHS stratifies each dimension individually. In this regard, 
the proposed Partially Stratified Sampling (PSS) methodology represents the broader class of all 
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possible stratifications between (and including) these extremes in which stratification occurs on 
low dimensional orthogonal subspaces of the A^-dimensional sample space S. As such, it can be 
viewed from either perspective as a compromise on true SS or as a generalization of LHS. Using 
either interpretation, true SS and LHS are special cases of PSS such that stratification occurs on N- 
dimensional and one dimensional spaces respectively. The following presents the PSS methodology 
and its properties in detail. 


3.1. PSS Method 

Let @i, i = 1,..., Ng denote Ng disjoint Aij-dimensional orthogonal subspaces of the A^-dimensional 
sample space S {Ni < N Vf) such that = S and &p n @q = 0 ',p ^ q. Note that the general 

form does not require Ni = Nj-,i ^ j. PSS divides each subspace 0* into a collection of Mi disjoint 
subsets Ptik'., /c = 1, 2,..., Mi. Lower Aj-dimensional random samples = {xiki,Xik 2 , ■ ■ ■, XikNi} 
are generated within each stratum Clik of subspace 0j according to the stratified sampling method 


in Section 2.2 Full N-dimensional samples x are assembled by randomly grouping the lower¬ 
dimensional samples generated in each subspace. That is, a low-dimensional sample Xj^ is randomly 
selected from each subspace 0* and these terms are grouped to produce a sample. 

Consider the simple case of a 4-dimensional sample space. Possible PSS subspace designs are 
shown in FigureNext, consider drawing four samples from the 4-dimensional sample space using 
a PSS 2x2 as shown in Figure First, we draw four 2-dimensional stratified samples in each 
subspace 0i and 02 (Figure]^). The 2-dimensional stratified samples are then randomly paired 
as shown in Figure [^. One could generalize this to cases where the strata are not equally sized in 
each subdomain but this will not be considered here given the challenges associated with computing 
sample weights. In particular, this represents the generalization of the Latin hyperrectangle sam¬ 
pling method of Mease and Bingham |25j in which it is necessary to violate the desirable property 
that the sum of the probability weights equals one. 

The notation used to describe a PSS design identifies, in decreasing dimensional order, the 
definitions of the subspaces as shown in Figure For low dimensional problems, this can be 
characterized by simply stating the dimension of all subspaces (e.g. PSS—2 x 2). However, for high 
dimensional problems, it is convenient to identify the number Ki of Aj-dimensional subspaces using 
a superscript. In general, this is denoted as PSS—N^^N^^ ... N^^^ where Ni > N 2 > ■ ■ ■ > N]\fg. 
For example, a 100-dimensional problem with 20 4-dimensional subspaces and 5 2-dimensional 
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Figure 2: Possible PSS subspaces for a 4D sample space. 


subspaces is denoted by PSS—Note that this alone does not fully specify a PSS. Additionally, 
it is necessary to stipulate which variables are being grouped in each subspace along with the number 
of samples and the number of strata in each subspace. 


3.2. Bias and Response Variance 

The response variance from PSS is derived here and shown to follow a nearly identical form to 
that of LHS. The derivation follows directly from [T] as given below. 

Partially Stratified Sampling divides an A^-dimensional sample space into Ng disjoint sub¬ 
spaces 0j. Each subspace is stratified into regions possessing probability volume —. Consid- 

n 

ering all possible combinations of strata spanning the subspaces yields Af-dimensional cells 
each having probability volume Each cell can be labeled by a set of Ng cell coordinates 

rrij = {mii,mi 2 ,... where rriij is the stratum index of cell i in subspace ©j. A Partially 
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Figure 3: Four samples drawn from a 4D sample space using a PSS 2x2. (a.) Draw stratified samples within the 2D 
subspaces, (b.) Randomly pair the 2D samples to assemble 4D samples. 


stratified sample of size n is obtained by randomly selecting n of the cells mj;i = 1,... ,n such 
that for each subspace &j, the set is a random permutation of the integers 1,... ,n. A 

random sample is generated in each selected cell. 

The probability density function of X given X G cell i equals f{x) if x G cell i and zero 
otherwise. The marginal distribution of Yi generated from a PSS is given by: 


p[Yi <y] = Y. ^ y\'^ ^ ^ 

cells 

9=1 


E /cell, 

f{x)dx 


-Ns 


cells h{x.)<y 

9=1 


lh{x)<y 


( 7 ) 


which is clearly equal to the distribution of Yi generated from SRS. Thus, a statistical estimator, 


































Tp, constructed from a PSS is an unbiased estimator of r. 

To derive the form of the variance of Tp, dehne the indicator variable a, as: 


di = 


1, if cell i is in the sample 
0, otherwise 


( 8 ) 


The PSS estimator can be written as: 


^ n ^ 


i=l 


where Yi = h(X.i) and Xj G cell i. The variance of Tp is given by: 


(9) 


Var[rp] = ^ ^^Ya.T[aig{Yi)] + ^ Cov[aig{Yi),ajg{Yj)] 


i=l 


i=l j=l 


Expanding the hrst term of Eq. (10) yields: 


( 10 ) 


YaviaigiYi)] = E[a^]Yar[giYi)] + E[giYi)fYar[ai] 


( 11 ) 


Exploiting the fact that E[ai] = E[af] = yields: 

s J- 


n" 


n ^ ^ n ^ 

J2YaT[aig{Y)] = :^JE^'^E[g{Yi) - + 


2=1 


2=1 


1 


1 


nY-i n2Ns-2 J 

2=1 




where gi = T[gf(y)|X G cell i\. By definition: 


( 12 ) 


E[g{yi) - yi? = / {g{y) - T?f{x)dx + {m - t? 

^cell 2 


which reduces Eq. (12) to: 


^ / 1 1 \ 
Yiiv[aig{Yi)] = nVar[y] - ^(//^ - r)^ + f j ^ 

i=l i=l ^ ^ i=l 


Consider now the covariance term in Eq. (10), which can be expanded as: 


(13) 


(14) 


E E Cov[aig{Yi), ajg{Yj)] = X] E yi^jEldiOj] - Yl yi^j 


i=l j=l 


i=l j=l 


i=i 

jy* 


(15) 
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The value of E[aiaj\ depends on the cell coordinates nij. If cell i and cell j have no cell coordinates 
in common then: 


E\aiaj\ = E[aiaj\aj = 0]P[aj = 0] + E[aiaj\aj = l]P[aj = 1] 

1 ( 16 ) 

(n(n — 

The case where cell i and cell j share at least one common cell coordinate corresponds to an 
inadmissible case (i.e. ai and aj cannot both equal 1) and thus: 


E[aiaj\ = 0 . 


(17) 


Combining Eqs. (14) and (15) and utilizing the relations in Eqs. (16) and (0 gives the following 
relation for the variance of a statistical estimator using PSS: 


Var[7p] = ^Var[T] - - rf + ( 

i=l ^ 


1 


1 


ry^Ns 




2=1 


+ 


n 


Ns-l 


(n- 1) 


Ns- 


R i=l j=l 


where R denotes the space of n^‘‘{n — 1)^“ admissible cell pairs with cell mean values 
Finally, utilizing and simple algebra, Eq. (18) reduces to: 

n — 1 1 

R 


VarfTp] = Var[rp] + 


n 71^“ (n — 1) 


-t). 


(19) 


As previously stated, the variance derived above follows directly from [T] with only slight mod¬ 
ifications based on the use of Ng stratified subspaces rather than N stratified lines. It follows 
precisely the same form and is straightforward to show that, when Ng = N, each stratified sub¬ 
space corresponds to a component of X resulting in a Latin hypercube sample. Thus, as in LHS, 
the proposed PSS method produces a variance reduction over SRS only when: 


n — 1 1 

n n^‘>{n — l)-^‘> 


R 


t) < 0. 


( 20 ) 


In other words, a variance reduction is achieved when the covariance between cells having no 
common cell coordinates is negative. 
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3.3. Asymptotic Properties 

Stein [3] provides proofs of several important asymptotic properties of Latin hypercube samples. 
Most notably, he shows that, in the limit n —)• oo, statistical estimates from Latin hypercube samples 
are normal possessing variance that is at least as small as an estimate using SRS. Furthermore, he 
proves that Latin hypercube sampling has the effect of hltering the additive components of h(X). 
In this section, we generalize these asymptotic properties to partially stratified samples. 

Consider the general transformation /i(X). By writing the transformation in the form: 

h(X) = h,(X)+r(X) (21) 


where /ia(X) is the best additive fit to /i(X) having form: 

N 

hai'^) = '^gi{xi) ( 22 ) 

i=l 

and r(X) is a non-additive remainder function describing variable interactions, Stein [3] shows that 
Latin hypercube sampling filters out the contribution ha(X.). Thus, the primary contributor to the 
variance of the estimate is the interaction function ^(X). 

We submit the following alternative form for /i(X): 

MX) = /(X) + r(X) (23) 


where /(X) takes instead the form: 

Ns 

/(X) = j;/,(X,) (24) 

i=l 

such that UjXj = X and rijXj = 0. In other words, instead of decomposing /i(X) into an additive 
sum of functions on a single variable, we decompose it into an additive sum of multivariate func¬ 
tionals. These multivariate functionals may or may not be additive - or even close to additive. If 
functional /i(Xj) is nearly additive, then it may be approximated as: 

N, 

/,(X,)«^5fc(xfc) (25) 

k=l 


and Latin hypercube sampling of this subspace is appropriate to filter out the contributions of 


9k{xk). When all such /i(Xj) Vi are closely approximated by the form in Eq. (25) then Latin hy¬ 
percube sampling of the full space will be effective. However, if any /j(Xj) is not well-approximated 
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by Eq. (25), alternative sampling methods should be considered on that subspace. In particular, if 


there are strong interaction terms it will be demonstrated that a stratified design of the subspace 
is more effective than a Latin hypercube for low-dimensional subspaces. 

An alternative presentation that leads to the same conclusion follows. Consider an iV-dimensional 
PSS possessing Ng orthogonal subspaces Qi]i = l,...,Ns each having dimension iVj such that 
For each subspace, define a new variable as follows: 


= /i(X,) 


(26) 


where X* is the vector of variables spanning subspace 0*, such that we can write Y = /i(X) ~ E(Z). 
In general, the functional forms of /*(•) and F(-) may be difficult or even impossible to determine. 
But, defining these variables conveniently reduces each subspace to a single dimension; effectively 
transforming a PSS on X into a Latin hypercube sample on the lower Xs-dimensional vector Z. 
Next, consider two cases: additive Zi and non-additive Zi. If Z, is additive then: 

Ni 

Zi = MXi) = Y,MXk) (27) 

k=l 

In this case, the most sensible approach is to perform Latin hypercube sampling on Xj - which is 
the same as defining Ni 1 -dimensional subspaces as opposed to a single Aj-dimensional subspace. 
In other words, higher-dimensional subspaces should not be applied to these variables. 

If, on the other hand, Zj is non-additive then performing LHS on Xj is not necessarily the best 
approach. Consider the following non-additive form for the variables Zj: 

Ni 

Z^ = lllkXt (28) 

k=l 

where the subscript k refers to each dimension of the subspace and 7 ^ and f3k are defined arbitrarily. 
Clearly, Zj are non-additive in but, it should be emphasized that this is not a general form for 
non-additive functions. It is used purely for convenience of demonstration. As will be seen, a true 
stratified design provides superior variance reduction on variable Zj when the dimension of Xj is 
low. 

Considering that a PSS can be essentially reduced to a Latin hypercube on the variables Zj, it 
follows that the asymptotic properties of partially stratified sampling can be inferred directly from 
those of Latin hypercube sampling given by Stein [3]. First, asymptotic normality is preserved 
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{Z — i are independent) - as is the property that a PSS will not increase the variance of a statistical 
estimator over SRS given a large enough N. Moreover, PSS has the effect of filtering out additive 
combinations of non-additive functions. In other words, a partially stratified design filters out the 
additive components in Zj. This is an important result because it states that, for a general non¬ 
additive function that can be written as an additive combination of simpler non-additive functions, 
the variance can be further reduced from a LHS by improving the sampling method specifically 
on the non-additive components. This will be demonstrated to have significant effect on variance 
reduction in the following sections. 

4. When to use partially stratified sampling? 

In this section, we explore the use of PSS. First, we discuss the benefits of true stratified 
sampling for low dimensional problems with strong interactions - showing that it produces a variance 
reduction that is superior to LHS and indicating that SS on low-dimensional subspaces where 
variables interact is likely to be preferable. Next, we explore the relationship between interaction 
strength and its influence on the decision to perform true stratified sampling or LHS on a subdomain. 
Lastly, we identify and discuss the primary challenges associated with selecting and implementing 
a PSS design. 

4-.1. Benefits of true stratified sampling in low dimension 

True stratified sampling has the benefit of improved variance reduction for low dimensional 
transformations possessing interaction terms as recently demonstrated by Shields et al. |22j . As 
illustrated in Eqs. ([^ and Q, SS and LHS reduce variance through different statistical mechanisms. 
SS reduces variance by considering the global statistical estimate to be the statistical mean of a 
set of local means conditioned upon the strata. In this way, variance is reduced for the main and 
interactive effects in equal measure. LHS, on the other hand, reduces variance by creating negative 
covariance between the transformed samples. As Stein shows, the strength of the covariance directly 
relates to the strength of the main effects (see [3], Theorem I). Thus, the stronger the main 
effects, the greater the variance reduction. However, for transformations possessing interactions, 
the variance reduction may be less than that achieved by SS. This is demonstrated clearly by the 
following simple example. 
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(a) (b) 

Figure 4: Variance of Monte Carlo estimates of the mean value for (a) an additive function and (b) a product function 
using stratified sampling (SS), Latin hypercube sampling (LHS) and simple random sampling (SRS). 

Consider two simple transformations. The first is an additive function defined by: 

2 ^ 

U = (29) 

i=l 

where Xi ~ 17(0,1) such that E\Yi] = 1. The second is a product function dehned by: 

N 

Y2 = X{X, (30) 

where Xi U (—\/3, \/3) such that E\Y 2 ] = 0. Note that, when E[Xi] / 0, Y 2 will also possess 
noticeable main effects although, as shown in [22], stratified sampling is still shown to reduce 
variance more than LHS in this case. In both functions, Xi are iid so it follows that Yi and Y 2 
are asymptotically normal and lognormal respectively as Ai —>■ 00 . Figure shows the variance of 
estimates for E\Yi] and E\Y 2 ] from 1,024 samples generated using SRS, LHS, and SS for different 
values of N. For Yi, LHS produces by far the best estimates of /iy^ (as expected based on Stein’s 
findings) yielding variance that is effectively negligible (~ 10“^*^). This is reduced by more than 
four orders of magnitude over SS, which yields still significant variance reduction over SRS. The 
variance reduction improves slowly for LHS with the dimension N but the variance of SS increases 
slowly with dimension given the inability to simultaneously stratify many dimensions. 

The product function, Y 2 , however tells a different story. For low-dimensional problems, SS 
produces a significant reduction in variance compared to LHS (which, as expected from Stein |3| 
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yields no variance reduction at all) that is up to three orders of magnitude. The reduction in variance 
over SRS/LHS diminishes with dimension although it maintains some minor advantage through 
moderate dimensions {N = 10). The performance loss results from the increase in variance of the 
strata means fj,i as the dimension grows. More specifically, although stratum probability pi remains 


constant, the extents of the strata in each direction grow as the dimension grows. For example, 

given 1,024 samples and strata probabilities pi = the stratum extents in a 2-dimensional 

2 


1024 


1 


space are ^ ^ = 


32 'V32 


1 


1 


1024 


= Pi) while the corresponding extents in a 10-dimensional space are 


1 Yl 


= 


10 


1 


= Pi)- 


2 J 1024 

These results make a compelling case for the use of PSS. For a general, high dimensional 
transformation possessing both main effects and interactions it will be advantageous to isolate the 
main effects and interaction effects to the largest extent possible and sample them in different ways 
(LHS for main effects and SS for interactions). The PSS methodology provides this capability. 


4-2. Effect of interaction strength 

The analysis above presents a difficult question: When is it appropriate to use true SS on a 
subspace and when should LHS be used? To investigate this question, we employ global or variance- 
based sensitivity analysis [26j to make a connection between interaction strength and the the choice 
of sampling method for a given subdomain. For this purpose, let us define the main effects of the 
transformation through the conditional variances: 


Fi = YaT{E[Y\Xi]} 

and the interaction effects through the conditional variances: 


(31) 


Vij =Yai{E[Y\Xi,Xj]}. 


(32) 


Using this notation, we can write the total variance VarjU} using the ANOVA-HDMR (Analysis 
of Variance - High Dimensional Model Representation) if X are independent variables as |26j : 


VT = J2v, + J2J2^^J + ■■■ + vu...k 

i i j>i 


(33) 


Such decomposition of variance is common in global sensitivity analysis where the sensitivity indices 
(Sobol indices) Si = Sij = etc. provide normalized contribution factors for the main effects 

Yrp Yj, 

and interaction effects respectively 
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Consider now the functional /i(Xj) from Eq. (24) and its decomposition into main effects and 


interactions as: 


/,(Xi) = /i,(X)+ri(X) 


(34) 


where /ja(X) and ri(X.) are the main and interactive contributions respectively. As a general rule, 
Latin hypercube sampling will produce significant variance reduction on /ja(X) and no variance 
reduction on rj(X). As shown by Stein [3], the variance of the LHS estimator Tn on this subdomain 
can be expressed as: 

Var[rLj] = “ y ri{x)^dFi{x) + 0{n~^) (35) 

which equates to standard Monte Carlo sampling on the interaction functional ri(X). Utilizing the 
ANOVA-HDMR decomposition, it follows that: 


n{xfdFi{x) = ^^Vij 
i j>i 


Vi2...k 


(36) 


yielding: 


Var[rii] 


X] X] -^ ^ 12 . ..fc 

i j>i 


+ 0(n-i) 


(37) 


such that the performance of LHS on a given subdomain can be directly related to the interaction 
sensitivities. 

Having established that true SS is more effective than LHS at reducing variance associated 
with interactions but is much less effective at reducing variance associated with main effects, we 
can make the qualitative statement that true SS is appropriate on a given subdomain when the 
contribution of the main effects to the variance from SS is smaller than the variance from LHS 
resulting from interactions. To elucidate this point, consider the simple example of a 2D quadratic 
functional with interactions: 


Zi = /i(Xi) = Xl + A| + CA 1 X 2 (38) 

such that it is reasonable to assume /ia(X) = and rj(X) = 0 X 1 X 2 when E'[Xi] = E\X 2 \ = 0 

(note the product X 1 X 2 may have considerable main effects when E[Xi]^ E[X 2 \ / 0). We consider 
two different distributions for X such that X X(0,1) and X ~ U{—V3,\/3) and study the 
effect of the coefficient c on the variance of the estimate of the mean E[Zi]. For the two 
distributions. Figure shows the interaction sensitivity indices computed using a Monte Carlo 
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Figure 5; Dependence of interaction sensitivity on the interaction coefficient, c, for a second order polynomial. See 
Eq. 


approach [26] as a function of c follows a logistic form when c is plotted in log scale. This suggests 
that interaction sensitivity is essentially invariant for small c (c < 0.1 corresponds to essentially 
negligible interaction) and large c (c > 10 corresponds to essentially negligible main effects) and 
varies strongly in the range 0.1 < c < 10. Figure further shows that /i(Xj) is more sensitive 
to interactions when X is uniformly distributed. Figure shows the variance of the estimate ^Zi 




Figure 6: Dependence of variance on the interaction sensitivity for stratified sampling and Latin hypercube sampling 
given (a.) normal input variables and (b.) uniform input variables. A very small degree of interaction dramatically 
increases the variance of LHS estimators relative to true SS. 


from LHS (variance essentially results entirely from interaction effects - see Eqs. (35) and (36)) and 
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true SS as a function of interaction sensitivity S 12 for the normal distribution (Figure |6(a)[ ) and 
uniform distribution (Figure [6(b)[ ). Also shown in Figure]^ is the variance of the stratified estimator 
attributed to main effects. These plots ratify our previous assertion that true SS is appropriate 
when the variance of the SS estimator associated with main effects is less than the variance of 
the LHS estimator associated with interactions effects but also highlights several other interesting 
and important factors in sample design selection. 1. For true SS, the variance associated with 
interaction effects is negligible unless the the interactions are very strong (S'12 > 0.5); 2. For LHS, 
the variance of the mean increases dramatically (orders of magnitude) with the introduction of 
even very small interaction effects. 3. SS becomes superior to LHS at very low sensitivity values 
(S'12 = 0.005 for uniform X and S12 = 0.1 for normal X). This leads to the conclusion that SS on 
a subdomain will likely be superior to LHS unless the interaction effects are weak. 


4-3. Challenges 

The most obvious and signihcant challenge to partially stratified sampling is to identify the 
optimal subspace decomposition. In general, the transformation /i(X) is rarely given in closed 
form and usually takes the form of a numerical computer model possessing complex interactions 
of variables. In certain cases, it may be clear which variables are interacting which will inform the 
PSS subspace definitions - but such cases are the exception. This challenge will be addressed in 
the following section. 

A second challenge is that, in order to be effective, the subspaces ©j must remain relatively low¬ 
dimensional. The true SS that occurs on these subspaces may perform poorly in high-dimensional 
spaces (such that e.g. 2^* > N). Consequently, it may be difficult to reduce the variances resulting 


from interactions of many variables as may occur if, for example, a component Zi in Eq. (28) is 
defined by the product of a large number of input variables Xj. 


5. “Latinized” Stratified Sampling 

It is possible to simultaneously reduce variance associated with the main effects and interaction 
effects by constructing, on a given Aj-dimensional subspace, a true SS that is at the same time 
a LHS. This is a achieved through a simple procedure referred to herein as “Latinized” Stratihed 
Sampling (LSS) (colloquially called “Su-Do-Ku Sampling” as a low dimensional LSS grid resembles 
the popular number game - see Figure]^. 
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Figure 7: Schematic description of the Latinized stratified sampling (“Su-Do-Ku Sampling”) method for 9 samples 
in 2-dimensions. 


The LSS method proceeds as follows and is shown graphically in Figure 

1. Draw a LHS from the subspace - do not group the individual variables (these are shown by 
an ‘x’ on the ordinates in Figure]^. 

2. Stratify the domain as desired ensuring that the stratification is consistent with an LHS design 
(solid lines in Figure]^. 
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3. For each stratum, randomly select a point Xi from each component of the LHS (without 
replacement) such that the sample x = {xi, X2, • • •, xn } lies within the stratum. 

4. Repeat for each stratum of the design. Notice that, after each sample grouping (i.e. after a 
sample is drawn from a given stratum), a number of LHS cells are now “off-limits” for further 
sampling. These cells are shaded in Figure 

5.1. Relation to Orthogonal Array-based LHS 

The LSS methodology is not an entirely new development. In fact, an LSS is exactly equivalent 
to an Orthogonal Array-based LHS (OA-LHS) [271 fT 6 ] when the strata all have equal probability. 
A subtle argument is that the LSS approach is perhaps more intuitive than the OA-LHS approach. 
Certainly, an LSS is easier to construct as it does not depend on generation of large orthogonal arrays 
that can pose significant challenges. Under the conditions of equivalence between LSS and OA-LHS, 
the LSS inherits the variance reduction properties that have been theoretically formulated by Tang 
m- That is, decomposing the transformation Y = L(X) according to its main effects {hi{Xi) = 
E[h{yL\Xj\ — fi) and bivariate interactions {hij{Xi, Xj) = E[h{X\Xi, X — j] — gi — hi{Xi) — hj{Xj)) 
as: 

N N i-1 

Y = ^hi(Xi) + h,j{Xi,Xj) + riX), (39) 

2=1 2=1 j = l 

the variance of an OA-LHS/LSS is given by: 

N N 

Var{ro} = Var{rR} - Var{/, (A,)} - - Var{^ (X„ X,)} + 0{n-^) (40) 

j=l i<j 

The novelty of the proposed approach is that the LSS can be applied to subspaces of different 
dimension in a simple and efficient way. Therefore, it is capable of targeting specific interactions and 
reducing variance associated with interactions of different order for different variable combinations 
without the need to generate large orthogonal arrays. For example, it is capable of reducing variance 
associated with main effects and interactions among, say, three variables in one subspace while 
simultaneously reducing main effects and interactions among two variables in another (disjoint) 
subspace. 

5.2. Challenges revisited 

The use of Latinized Stratified Sampling in combination with Partially Stratified Sampling 
(LPSS) mitigates the challenges posed in Section 4.3 to a large degree. The primary challenge with 
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PSS is deciding upon a set of appropriate subspaces on which to sample. Using PSS in its basic 
form, the decision to stratify two variables together based on an assumption that they interact will 
result in a significant variance cost if those two variables do not actually interact. Meanwhile, the 
opposite choice will have a similar cost if they do interact. With the introduction of the LPSS, 
it is sufficient to stratify a set of variables together simply based on the possibility that they may 
interact. If they do, in fact, interact then the savings will be amplified by reduction in both the 
main effects and the interactions. If they do not interact, there will be no increase in variance since 
the main effects are also being filtered. 

6. Demonstration Problems 


In this section, we present a series of high dimensional problems with significant interactions in 
order to demonstrate the proposed methodology and its potential benefits and shortcomings. 

6.1. High dimensional polynomial functions 

Consider a general high dimensional second order polynomial of the form: 




Ki 


Kj, 


y — ^ oikXl + ^ l3kX2k-iX2k + ^ 7kXk 


(41) 


k=l k=l k=l 

where Ki\f 2 refers to the number of non-interacting second-order terms, Kj refers to the number of 
second-order interaction terms, and Kj\fi refers to the number of non-interacting first-order terms. 


The dimension of the polynomial expressed in Eq. (41) is equal to K = m.ax{K]\f 2 ,‘^Hj,K]sri). 
In the following we consider problems of this form with dimension K = 100 and two different 
distributions: (a.) Xk ~ A^(0,1); and (b.) Xk ~ A^(l,l). The shifted mean has the effect of 


introducing main effects into the second (interaction) term of Eq. (41). 


Several cases of the polynomial model given in Eq. (41) are studied as outlined in Table 1. In 
each case, we consider ak = fdk = Jk = 1, V/c. Eor practical purposes the relative influence of 
each term, as specified by these coefficients, will play an important role in the response Y, and by 
extension will influence the decision to use PSS versus LHS and how to define the PSS subspaces. 
However, in the interest of conciseness and considering the short study presented in Section 4.2, we 
focus on the inffuence of the existence of the various terms and not on their relative strengths. 

Six different sampling methods have been employed for each case identified in Table 1. In each 
case 625 samples are drawn. As a baseline for comparison, both SRS and LHS are employed. Two 
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Case 

Kn2 


Ki 


Kni 

1 

100 

0, 1, 

2, 5, 10, 

25, 50 

100 

2 

100 

0, 1, 

2, 5, 10, 

25, 50 

0 

3 

0 

0, 1, 

2, 5, 10, 

25, 50 

100 

4 

0 

0, 1, 

2, 5, 10, 

25, 50 

0 


Table 1: Dimension and interaction terms for the high dimensional polynomial test problem given by Eq. ( |41[ ). 

different variations on PSS and LPSS are also used. The first utilizes 2D subspaces for all partial 
stratification (denoted PSS-2^° and LPSS-2^° respectively). The second utilizes 2D subspaces for 
interacting variables and ID subspaces for non-interacting variables (denoted PSS-2^^ ■^k- 2 Ki 

LPSS-2^n'^-2'^0- 






(d) Case 4 


Figure 8: Variance of Monte Carlo estimates of the mean value for a second-order polynomial function given in Eq. 
(|41[|. Input variables are N(0,1). Description of the cases are given in Table [TT] 
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Figure 9: Variance of Monte Carlo estimates of the mean value for a second-order polynomial function given in Eq. 
(411. Input variables are N(l,l)- Description of the cases are given in Table 6.1 


The variance of the E\Y] estimate from the 625 has been computed for each sampling method 
from 1000 Monte Carlo simulations. Figure plots these variance estimates for cases 1-4 as a 
function of the number of interaction terms Kj given Xk ~ A^(0,1). The same plots for Xk ~ 
A^(l, 1) are provided in Figure]^ From these plots, several interesting features are observed. 

• LPSS-2^^1^“^^^ and LPSS-2^^ are consistently the best methods in terms of variance reduc¬ 
tion. Both perform nearly identically and provide a very large variance reduction - especially 
when there are both strong main effects and interactions. 

• The variance reduction from PSS-2®*^ depends strongly on the relative strength of the interac¬ 
tions and the main effects. When the main effects are significant PSS-2®*^ is often less effective 
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than LHS. The reason for this is that it stratifies the space based on interactions that are 
not present at the expense of the main effects in those variables. However, as the interaction 
strength grows {Kj gets larger), PSS-2^° begins to reduce variance compared to LHS. The 
point at which this occurs depends on the strength of the main effects. For example, in Cases 
3 and 4, where the main effects are weak, PSS-2^^ is nearly always superior to LHS. On 
the other hand, in Cases 1 and 2 where the main effects are strong, PSS-2®*^ only becomes 
advantageous for a large number of interaction terms {Kj = 40 — 50). 


As expected from Stein [3], when no main effects are present, LHS does not reduce variance 
(N(0,1) Case 4 - Fig. 8 (d)| ). However, all PSS and LPSS methods provide considerable variance 
reduction in this case. 


6.2. Rosenbrock function 

The second problem we consider is the Rosenbrock function [28] defined as follows: 

K-l 

•Pr(x) = ^ [l00(x- - Xi+if - {xi - 1)2] (42) 

i=l 

where x = {xi,X 2 ,... ^xk} is a AT-dimensional random vector. The Rosenbrock function is a 
strongly nonlinear high-dimensional function that is commonly used as a benchmark problem for 
high-dimensional optimization. We consider the case where K = 100 and Xi 17(0,1). 

Noticing that, unlike the polynomial function defined in the previous section, although it has 
only bivariate interactions, these interactions occur between all subsequent variables (e.g. term 
1 interacts with term 2 while term 2 independently interactions with term 3). Thus, the choice 
of partial stratification is less obvious. We consider two different PSS/LPSS designs that group 
two terms (PSS—2®^) and four terms (PSS—4^^) respectively. The mean value E[Fji\ = 2013 is 
estimated from 625 samples and the standard deviation is computed by repeating the estimation 
of E[Fji\ 5000 times. The results using four different sampling methods (SRS, LHS, PSS, LPSS) 
are given in Table 2. 

Notice that, while LHS provides some variance reduction, the PSS and LPSS methods are 
considerably more effective. Notice also that the 4D partial stratification reduces variance more 
than the 2D partial stratification. This is because the higher degree of stratification is able to reduce 
variance associated with more interactions. Finally, similar studies conducted over different ranges 
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Method 

Standard deviation of mean estimate 

Rosenbrock 
E[Fr] = 2013 

Schwefel—X(0,1) 
E[Fs] = 5057 

Schwefel—X(l, 1) 
E[Fs] = 343,400 

SRS 

8.778 

247.5 

3131.0 

LHS 

6.756 

234.1 

245.9 

PSS-2'50 

4.856 

227.1 

358.5 

PSS-425 

4.588 

234.2 

983.7 

LPSS-250 

4.819 

220.1 

241.0 

LPSS-425 

3.813 

226.9 

236.2 


Table 2: Standard deviation of mean value estimates for the Rosenbrock and Schwefel’s functions. 


of the Rosenbrock function (Xj ~ U{1, 2), Xj ~ C/(0, 2), Xj ~ C/(0, 3), Xj ~ C/(0, 5), Xj ~ [7(9,10)) 
yield trends consistent with those in Table 2. 

6.3. Schwefel’s Problem 1.2 

The final numerical demonstration problem is the Schwefel’s problem 1.2 that possesses very 
high dimensional interactions defined by: 

K 


Fs{^) = Yl Y.'- 


(43) 


\i=i 

Again, we consider K = 100 and consider the standard deviation of estimates of the expected 
value from 5000 repeats of 625 samples. As summarized in Table 2, we consider Xj ~ X(0,1) and 
Xj ~ X(l, 1) with true mean values estimated from Monte Carlo simulation of 5,057 and 343,400 
respectively and two different partial stratifications (PSS/LPSS—2®° and PSS/LPSS—4^^). In the 
case of Xj ~ X(0,1), it is clear from Eq. ( |43| ) that the main effects are very small. Hence, not only 
is the standard deviation very large (yielding a coefficient of variation on the estimate of nearly 
5%), but LHS produces essentially negligible variance reduction. Similarly, the PSS and LPSS 
methods are equally ineffective at reducing variance in this problem. However, given a shifted 
mean Xj ~ X(l,l), expansion of Eq. (43) yields significant main effects. Thus, LHS produces 
more than an order of magnitude variance reduction. PSS is slightly less effective since it is able 
to reduce variance associated with the <2nd and <4th order interactions respectively but cannot 
reduce variance associated with the significant main effects. The LPSS method remains effective 
because it reduces variance associated with both effects. 

This problem serves to illuminate the limitations of both the proposed methods and existing 
variance reduction techniques. There is currently no effective means of reducing variance associated 
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with very high-dimensional interactions. Fortunately, such high-dimensional interactions typically 
occur only in select circumstances. 


7. Application to plate buckling 


The assessment of structural component strength often involves the interaction of various ma¬ 
terial and geometry parameters. One such case is the analysis of plate buckling strength which is 
of interest in many structural applications. Consider a rectangular plate that is simply supported 
on all four edges subjected to uniaxial compression. Analytical expressions have been derived to 
assess the buckling strength, first by Faulkner |29j, who derived the normalized buckling strength 
for a pristine plate as: 


(To V A2 

where is the ultimate load capacity, (Tq is the yield stress of the material, and 


(44) 


. _b [o^ 
~ t\ E 


(45) 


is referred to as the slenderness of the plate with width b, thickness t, and elastic modulus E. 
This equation was later modified by Carlsen m to include the effects of non-dimensional initial 
deflections 6q = ^ where wq is the magnitude of the deflection and residual stresses resulting from 
welding at the perimeter as: 


(/> = 





(46) 


where rjt is the width of the zone of tension residual stresses. 

Uncertainty analysis of such plates has long been of interest to the naval community m and it is 


clear from Eq. (46) that the various parameters controlling the buckling strength and its variability 
interact strongly. We are interested in defining the “best” stratification of the domain such that 
variance of response statistics from a Monte Carlo analysis of plate buckling will be minimized. 
We consider mild structural steel plates commonly used in naval applications with material and 
geometric variabilities estimated from the data presented by Hess et al. |32] and Guedes Soares m 
as presented in Table 3. The plate width b is estimated based on statistical analysis of variability 
in stiffener spacing for common ship plates. Thickness variations are estimated from statistical 
analysis of plate thickness from over 2000 measurements [32] while elastic modulus and yield stress 
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Variable 

Physical Meaning 

Nominal Value 

Mean 

COV 

Distribution Type 

b* 

width (in.) 

24 

0.992 X 24 

0.028 

Normal 

t* 

thickness (in.) 

0.5 

1.05 X 0.5 

0.044 

LogNormal 


yield stress (ksi) 

34 

1.3 X 34 

0.1235 

LogNormal 

E* 

elastic modulus (ksi) 

29,000 

0.987 X 29,000 

0.076 

Normal 

X '^0 

& = T 

7] 

initial deflection 

0.35 

1.0 X 0.35 

0.05 

Normal 

residual stress 

5.25 

1.0 X 5.25 

0.07 

Normal 


Table 3: Distributions for plate material, geometry, and imperfection variables. Distributions based on guidance from 
|32 | and |31| . 


result from statistical analysis of common structural steel such that the large mean value and 
COV for yield stress is a result of the common practice of downgrading higher strength steels that 
do not meet specification and the fact that Navy standards require only a minimum yield stress. 
Lastly, the initial deflections and residual stress variations are prescribed according to the findings 
of Faulkner |2^ and Antonia |3^ as reported by Guedes Soares m- 

As illustrated in Table 3, the plate buckling problem considered here is 6-dimensional. Thus, 
there are many possible partial stratifications. We consider three different PSS/LPSS designs (4 if 
LHS is included) as shown in Table 4 that attempt to group variables according to the interactions 


suggested by Eqs. (45) and (46). In each case, a total of 625 samples are drawn. For this problem. 


the variance reduction associated with hltering the main effects is more significant than that of 
filtering just the interactions. That is, LHS provides a greater variance reduction than all three 
PSS designs. However, the problem is able to benefit from reducing the variance associated with 
both main effects and interactions as all three LPSS designs provide a greater variance reduction 
than LHS, with LPSS—4^1^ providing the largest variance reduction. Thus, we conclude that the 
LPSS-—4^1^ is the most effective sample design for this problem. Referring to the discussion of 
Section 3.3, this design has a similar effect (though not exactly) to constructing a Latin hypercube 


sample on the three variables A, dp, and rj. Given, the form of Eq. (46), this is somewhat intuitive. 


8. Conclusions 


In this work, the Latin hypercube sampling method of [T] has been generalized by defining a 
spectrum of stratified sampling methods of which true stratified sampling and Latin hypercube 
sampling lie at its extremes. The intermediate designs on the spectrum are defined as Partially 
Stratified Sample designs and are shown to reduce variance associated with variable interactions 
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Design 

Paired Variables 

Strata 

Mean Strength E[(j)\ 

Std. Dev. of E[(j)\ 

SRS 

N/A 

N/A 

0.5590 

2.00e - 3 

LHS 

[b]- [t]; [aoj; [E]-, [5o]; [v] 

[625,625,625,625,625,625] 

0.5590 

4.33e - 4 

PSS-23 

[b,t]- [ao,E]; [(Jq,??] 

[25,25,25] 

0.5590 

6.36e-4 

PSS-2212 

[b,t]- [ctoW]; N; W 

[25,25,625,625] 

0.5590 

6.17e-4 

PSS-4H2 

[b,t,ao,E]; [^o]; W 

[5,625,625,625] 

0.5590 

l.lOe - 3 

LPSS-23 

[b,t]- [(7 o,E]; [(io,??] 

[25,25,25] 

0.5590 

3.87e - 4 

LPSS-2212 

[b,t]- [ao,E]- [5o]; W 

[25,25,625,625] 

0.5590 

3.95e-4 

LPSS-4H2 

[b,t,ao,E]; [5o]; [rj] 

[5,625,625,625] 

0.5590 

3.45e - 4 


Table 4: Distributions for plate material, geometry, and imperfection variables. Distributions based on guidance from 

[32] and |3Tj . 


whereas LHS reduces the variance associated with the main effects of the variables. The properties 
of PSS designs are studied and some considerations and challenges associated with their use are 
discussed. The challenges and shortcomings of PSS designs are mitigated by coupling a PSS design 
with a newly introduced stratified sampling method called Latinized stratified sampling that is 
simultaneously a LHS design and a fully stratified design. The LSS method is equivalent, under 
certain conditions, to an Orthogonal Array based Latin hypercube but is signficantly simpler to 
obtain. Merging the LSS method with the PSS method yields the Latinized Partially Stratified 
Samping method (LPSS) which is shown, for many high dimensional applications, to provide su¬ 
perior variance reduction when both low-order interactions and main effects are present. Several 
numerical examples have been provided and the new methods have been applied to a plate buckling 
problem in structural mechanics. 
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